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SUMMARY 


Second- and third-order two time-level five-point explicit upwind-difference 
schemes are described for the numerical solution of hyperbolic systems of conservation 
laws and are applied to the Euler equations of inviscid gas dynamics. Nonlinear 
smoothing techniques are used to make the schemes total variation diminishing. In 
this method both hyperbolicity and conservation properties of the hyperbolic conserva- 
tion laws are combined in a natural way by introducing a normalized Jacobian matrix 
of the hyperbolic system. Entropy satisfying shock transition operators (which are 
consistent with the upwind differencing) are locally introduced when transonic shock 
transition is detected. Schemes thus constructed are suitable for shock-capturing 
calculations. The stability and the global order of accuracy of the proposed schemes 
are examined. Numerical experiments for the inviscid Burgers' equation and the com- 
pressible Euler equations in one- and two-space dimensions involving various situa- 
tions of aerodynamic interest are included and compared. Numerical results using 
second-order schemes indicate that the stationary shock can be captured within (at 
most) two transition zones and the contact surface can be accurately resolved, and are 
comparable with those obtained by van Leer's and Harten's high resolution schemes. 
Preliminary results for one-dimensional problems using the third-order scheme are 
presented . 

1 . INTRODUCTION 


The hyperbolic conservation equations of inviscid gas dynamics have two mathe- 
matical properties which are of great significance in the design of a numerical 
method. First, for nonlinear hyperbolic partial differential equations, the theory 
of characteristics provides the correct directional signal propagation information and 
physical domain of dependence and has long been recognized as a natural procedure for 
solving the equations (refs. 1 and 2). Second, the conservation- law form of the equa- 
tions permits one to construct conservative shock-capturing numerical schemes 
(refs. 3-5) due to a theorem of Lax and Wendroff (ref. 3). 

A first-order accurate upwind finite difference scheme was developed for solving 
the nonlinear hyperbolic equations by Courant, Isaacson and Rees in 1952 (ref. 1). 
Their method was based on the normal or characteristic form of the quasi-linear first- 
order hyperbolic system. With the idea in mind of preserving the domain of depen- 
dence, the spatial derivative in each direction of the characteristic form was approx- 
imated by either a first-order forward or backward finite difference quotient depend- 
ing upon whether the local eigenvalue (characteristic speed) is negative or positive. 
The notion of eigenvalue splitting was introduced to provide the automatic switching 
for constructing uniform upwind finite difference schemes. 

Attempts to construct finite difference schemes by using both the conservation 
property and hyperbolicity property of the hyperbolic conservation equations were 
presented by Godunov (ref. 6) and by Steger and Warming (ref. 7). Recently, several 
upwind flux difference splitting methods have been proposed (refs. 8-12) for solving 



hyperbolic conservation equations. A review of some recent developments in one-sided 
(upwind) difference schemes has been reported by Harten, Lax, and van Leer (ref. 13). 

In references 14 and 15, a first-order upwind flux difference splitting method 
was developed for hyperbolic systems of conservation laws and was applied to the 
Euler equations of inviscid gas dynamics. First-order explicit and implicit schemes 
are described and can be viewed as the corresponding development for hyperbolic con- 
servation laws as the Courant-Isaacson-Rees (ref. 1) method for the hyperbolic equa- 
tions in quasi- linear form. The scheme has similar desirable shock calculation 
properties of Osher's scheme (ref. 10) but the integration path used in his method is 
avoided. In comparison, the special averaging technique used in Roe's scheme 
(ref. 8) was replaced by a simpler formulation. 

It is well known that first-order accurate schemes are too diffusive, but classi- 
cal higher-order schemes, while giving better resolution of the discontinuities of 
the solution, exhibit spurious oscillations around such points (e.g., the Lax-Wendroff 
scheme (ref. 3), the Fromm scheme (ref. 16), and the Warming-Beam scheme (ref. 17). 

In recent years, effort has been put into obtaining second- or higher-order schemes 
which give high resolution while remaining total variation diminishing or nonincreas- 
ing (TVD), a notion introduced by Harten (ref. 18). 

The present study is a sequel to the work contained in references 14 and 15. 
Second- and third-order (in space and time) accurate explicit two time-level five- 
point upwind difference schemes are described for the numerical solution of hyperbolic 
systems of conservation laws and are applied to the compressible Euler equations of 
inviscid gas dynamics. Nonlinear smoothing techniques (or flux limiters) are used to 
make the schemes "monotonicity-preserving" (or TVD). Many such nonlinear smoothing 
techniques are available. Recently, a systematic study on flux limiters including 
those proposed by van Leer (ref. 19), Roe (ref. 20), Harten (ref. 18), and 
Chakravarthy and Osher (ref. 21) has been reported by Sweby (ref. 22). Implicit TVD 
schemes have been presented by Harten (ref. 23) and Yee, Warming, and Harten 
(ref. 24). 

Section 2 begins with a brief review of the basic notion and the first-order 
characteristic flux difference splitting method described in references 14 and 15. In 
Section 3, shock transition operators which allow entropy-satisfying shock transitions 
are described in accord with the principle of upwind differencing. In Section 4, 
five-point second-order upwind difference schemes for single conservation laws and 
for systems are given. They are constructed based on Fromm's zero-average phase- 
error scheme (ref. 16). A nonlinear smoothing technique devised by van Leer (ref. 19) 
is used to make the scheme monotonicity-preserving. We also describe a second-order 
scheme obtained by applying our first order accurate scheme to a modified flux func- 
tion similar to the one devised by Harten (ref. 18). In Section 5, we extend the 
second-order upwind scheme using van Leer's (ref. 19) flux limiter to a third-order 
upwind scheme based on the characteristic flux difference splitting concept. Stabil- 
ity and accuracy of the second- and third-order schemes are examined in Section 6. 
Applications to the Euler equations of gas dynamics are illustrated in Section 7. In 
Section 8, numerical experiments for the inviscid Bergers' equation and the Euler 
equations in one- and two-space dimensions involving various situations of aerodynamic 
interest are considered and the results are discussed. Concluding remarks are made 
in Section 9. 


The author wishes to thank Bob Warming for many stimulating discussions and Helen 
Yee for bringing to his attention Sweby 's article. I am indebted also to Dennis 
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Jespersen for the valuable criticisms of the early version of this manuscript. The 
constant interest and encouragement of Harvard Lomax during the course of this work is 
gratefully acknowledged. This work was carried out under the auspices of the 
National Research Council. 


2. CHARACTERISTIC FLUX DIFFERENCE SPLITTING 


To illustrate the basic notions, we consider numerical solutions of a hyperbolic 
system of conservation laws in one space dimension 

3 U + 8 F = 0 (2-la) 

t x 

subject to the initial data 


U(x,0) = (|>(x) > -co < x < oo (2-lb) 

where <j> (x) is a given function and U = (u x ,u 2 , . . . , 1 ^)^ and F = (f 1 ,f 2 , . . . ,f m ) T 
are m-component column vectors. Equation (2- la) can be rewritten as a quasi-linear 
system 


3 U + A3 U = 0 (2-2) 

t x 

where A is the Jacobian matrix 3F/3U. We require that the sytem of equa- 
tions (2-la,b) be hyperbolic in the sense that the Jacobian matrix A has real eigen- 
values and a complete set of eigenvectors. Let the eigenvalues of A be denoted by 
a 1 , a 2 , . . . ,a m with corresponding right eigenvectors r x ,r 2 , . . . ,r m . Defining Q -1 (U) 
to be the matrix whose jth column is r j (U) , it then follows that 

QAQ -1 = A (2-3) 

where A is a diagonal matrix with the eigenvalues a^ as its elements and the 
norms of Q and Q _1 are uniformly bounded. By using equation (2-3), equation (2-2) 
can be put in the following characteristic form: 

Q -1 3 U + AQ _1 3 U = 0 (2-4) 

t x 

For the purpose of analysis, we assume that the coefficient matrix A is "frozen," 
that is, constant. By virtue of equations (2-3) and (2-4) one can transform equa- 
tion (2-4) to the uncoupled system: 

3 fc v^ + a £ 9 x v £ =0, £=1,2, ...,m (2-5) 

by defining a new characteristic vector V = (v 1 ,v 2 , . . . ,v m )^ = Q 1 U. The eigenvalues 
are real and equation (2-5) is equivalent to a set of uncoupled scalar wave equations. 
For simplicity, the subscript £ will be dropped for the remainder of this section. 

Define a uniform computational mesh {xj,t n }, with mesh sizes Ax and At. The 
discrete representation of v(x,t) on the mesh is Vj and A = Ax/At is the mesh 
ratio . 
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A first-order upwind difference scheme for solving equation (2-5) was introduced 
by Courant, Isaacson, and Rees (ref. 1) 


n+i n 

v . = v . 

J J 


f j-(l/2) Vn ’ 

Xa< 

b + (i / 2 ) vI1 5 


a > 0 
a < 0 


( 2 - 6 ) 


We can combine the above scheme into one uniform scheme by splitting the eigenvalue 
a into the positive part a + and negative part a~. 


n+i 
v . 

J 


n , + . n 

v j - Aa V(vu v 


Aa 


V(1/2) V 


(2-7) 


where 


a 


+ 

a 


+ a 


( 2 - 8 ) 


and 


+ 


a 


(a + |a|)/2 = max(a,0) 


a = (a - j a I ) /2 = min(a,0) 


(2-9) 


with 


A . , , N v 

3- (l /2 ) 


n n 

V . - V . 

J J- 1 


Extending the above procedure to system equation (2-2), one obtains 

yii+1 _ y n , .+ 4 TT n . „n 

j j 


AA - / / sA. “ AA . . . A . . .U" 

l-(i/2) j-O/2) J+(i h) J + (l/2) 


( 2 - 10 ) 


with 

A = A + + A" (2-11) 

where 

A 1 = QA ± Q -1 (2-12) 

and 

A = A + + A" (2-13) 

where 

A~ = diag{a*} (2-14) 

In references 14 and 15, we presented a corresponding development for hyperbolic 
systems of conservation laws equation (2-1) by using the following approach: For the 

purpose of upwind differencing, we will rewrite equation (2-la) as 

9 U + (A + + A") 9 F = 0 (2-15) 

t x 
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a + ^ i — x 

where A (U) = QA Q remain to be defined but are required to satisfy 

A + + A" = I (2-16) 

/s + ^ — _j_ __ 

A and A are closely related to matrices A and A of equation (2-11) and carry 
the hyperbolicity property of the equations through the expression 

A + = Qdiag{a^}Q _1 and A = Qdiag{ a^}Q -1 (2-17) 

The "normalized" eigenvalues are defined as 

± 

a l = = 2 ± \ sgn(a ^ * if a n * 0 (2-18) 

For a^ = 0, which corresponds to a sonic point, special treatment is needed. Suppose 
that the eigenvalue a^ is equal to zero at the nodal point J. Then in order to be 
consistent with upwind differencing we take 

aj = aj +1 if a j = 0 (2-19) 

This will exclude stationary expansion shocks and make the scheme entropy satisfying. 
Further discussion on this aspect will be given in the next section. 

An explicit first-order upwind difference scheme for equation (2-15) can be 
expressed as 


U n +l _ yn = a. F n - AA.. . .A.,. , s F n 

3 3 3"(i/2) 3-(l/z) 3 + (l/2) 3 + d/ 2 ) 

For steady-state calculations, a corresponding implicit scheme is given by 


( 2 - 20 ) 


[I + XA j-(l/2) A j-(l/2) Al1 + ^ +(l / 2 )V(l/0 An]AUn ^ 


( 2 - 21 ) 


of (2-20) , where AU n = U n+1 - U n . 


The nodal point at which the coefficient matrices are evaluated is relevant to 
the numerical schemes thus generated. We list two ways of defining the values of 
At. , / x used in the present study: 

J±(l/2) y 1 


(i)A + . , . = (A + 4- A + )/2 , A. , . , . = (A. + A. )/2 

3- ( 1 / 2 ) j-i 3 3 + (i/2) 3 3+i 


( 2 - 22 ) 


( ii ) AT 


d/2) = A ^j-i + V /2] ’ Vd/a) = A_ [ (U j + Vi )/2] (2 " 23) 


We wish to point out that the averages given by (2-22) or (2-23) are only correct for 
mesh intervals in which the eigenvalues are of the same sign. For mesh intervals 
which cover regions of transonic shock transitions (i.e., regions where eigenvalues 
switch signs) , the above formulas are not consistent with the idea of preserving the 
domain of dependence. 
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A + 

Examples for constructing matrices A for the Euler equations of gas dynamics 
will be given later. 

It is noted that the following relations hold for the normalized Jacobian 
matrix A~: 

A~A = A~ , and A“A 2 = A" (2-24) 

These relations can be used to reduce the computational effort in the second- and 
third-order schemes which are to be described later. For example, we can evaluate 
aA~ + bA'*’ + cA ±2 as Qdiag{aa^ + ba~ + caj }Q -1 . It is also noted that 

A ± 8 U = A~A3 U = A ± 9 F (2-25) 

XXX 

Symbolically, we can view "A~ = A~/A," which enables us to convert schemes for 
quasi-linear hyperbolic equations to schemes for hyperbolic conservation laws very 
directly and vice versa. 

Equations (2-10) and (2-20) are equivalent systems if no transonic shock transi- 
tions occur. One needs some information about the propagation speed when transonic 
shock transition takes place. Equation (2-10) does not have this mechanism, hence 
shock-fitting (or tracking) is needed. Equation (2-20) allows us to do shock- 
capturing calculations owing to its conservative property. It is observed that, in 
general, first-order schemes are adequate for the resolution of shock discontinuities; 
but for contact discontinuities, higher-order schemes are desirable. 


3. ENTROPY SATISFYING SHOCK TRANSITION OPERATORS 


The scheme defined by equations (2-20) and (2-22) is suitable for time-dependent 
calculations. It was observed that, for stationary shock calculation, spikes appear 
sometimes at the shock transition zones due to violations of the physical domain of 
dependence. It was also found that scheme (2-20) itself, when applied to the inviscid 
Bergers' equation, admits stationary expansion shocks (consider f(u) = u 2 /2 with 
initial data u 0 = -1, x< 0, u 0 = +1 , x > 0) and, therefore, is not entropy satisfying. 

The Engquist-Osher scheme has a mechanism to provide an entropy satisfying tran- 
sonic shock transition (spreading over at most two transition cells) through an inte- 
gration formulation (ref. 9). It was investigation of the scheme of Engquist and 
Osher which prompted the work of the present section. 

It is important to understand what kind of transition path a conservative, shock- 
capturing, finite difference scheme takes when differencing across a transonic transi- 
tion zone where the eigenvalue switches sign. Differencing across shock transition 
zones is necessary for satisfying the Rankine-Hugoniot jump relation and for generat- 
ing the correct shock speed, hence allowing one to do shock-capturing calculations. 
Most schemes do not have a special design to handle this transonic shock transition 
and this transition spreads over several "so-called" transition cells; sometimes this 
may damage the scheme. 
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In this section, we address this transonic shock transition problem for a conser- 
vative upwind difference scheme according to the upwind differencing principle and 
the entropy conditions. 

A shock transition operator is introduced into scheme (2-20) when transonic shock 
transition is detected. There are two types of transonic shock transitions: 

(1) transonic shock (where the eigenvalue or characteristic speed changes from posi- 
tive to negative) and (2) transonic expansion (where the characteristic speed changes 
from negative to positive) (see fig. 1). 

For illustrative purposes, we consider a single conservation law: 

3 u + 3 f (u) = 0 , a(u) = 3f/3u (3-1) 

L X 

Equation (2-20) for (3-1) becomes: 


. n n+i n , - + , r n ,-n , 

Au . = u. - u. = -Aa. . , .(f, - f. ) 
3 3 3 3-(i/2) 3 3-1 


„n 


Aa , / v (f . . 
:+(i/2) 3+1 




(3-2) 


The shock transition operators are applied when shock transition is detected. For 
example, if the shock transition lies between nodal points J and J + 1 (as depicted 
in fig. 1) , i .e . , 

a J a J+i < 0 (3-3) 

We can approximately locate the position of the sonic point, x* = (J + 0)Ax where 
a* = 0, by the following formula: 


6 = ” a J /(a J+i " a J> (3 " 4 > 

The state variable u* and the flux vector f* at the sonic point can be approxi- 
mated by 

U * = u j + e < u j +1 " u j> 

and 


f* = fj + ©(fj^ - fj) or f* = f(u*) 

After locating the sonic point, we can separate the interval [J Ax, (J + l)Ax] into 
two parts [J Ax,(J +0)Ax] and [ (J + 0)Ax,(J + l)Ax] by the sonic point. In each part 
the eigenvalue distribution is of the same sign, either positive or negative. Then 
we can apply upwind differencing according to the sign of the eigenvalue in each 
region which the flux difference interval covers. 

To update the state variable at nodal point J, we use 


. n n+i n , -+ c N 

Au T = u T - u = -Aa T . , . ( f T - f T ) 
J J J J-(l/2) J J-i' 


(a + a,) 


(f* " fj) ~ X 


( a* + aj +1 ) 


^ f J+i f *^ 


(3-5a) 
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and to update the state variable at nodal point J + 1, we use 


. n n+i 

Au t , = u T . 

J+i J+i 


U J+1 Aa j+(3/2) ^ f J+2 f J+l^ 


-X 


(a J + 


at) 


(f* - 


f J> ' 


/'+ , ^+ 
< a * + a 


J+i- 


(f 


J+i 


f*) 


(3-5b) 


In equations (3-5a,b), the split normalized eigenvalues a” need some special treat- 
ment in order to be consistent with the upwind differencing. Since the sign function 
for zero is indeterminate, one way is to use the asymptotic values (0 + ) and (0~). 

For equation (2-18), we have 

sgn(0 + ) = 1 and sgn(0 ) = -1 (3-6) 

To determine the values of a~, we take 


and 


+ + 



(3-7 a) 


(3— 7b) 


In so doing, we can avoid entropy-violating transition and the differencing is consis- 
tent with the principle of upwind differencing of Courant, Isaacson, and Rees 
(ref. 1). 


For a single conservation law, the scheme defined by equations (3-2), (3-5a), 
and (3-5b) is identical to the scheme of Engquist and Osher (ref. 9) which is also 
identical to the flux vector splitting scheme of Steger and Warming (ref. 7). The 
same procedure' can be applied to the systems (2-20). For example, in one-dimensional 
gas dynamics, the equation corresponding to equation (3-3) is 

(u - c ) j (u - c) J+i < 0 (3-8) 

where u is the fluid velocity and c is the local speed of sound of the fluid. 

Numerical examples for a single conservation law, and for systems of hyperbolic 
conservation laws involving moving shocks, stationary shocks, and transonic expansions 
are considered and the results are included in a later section. 


4. SECOND-ORDER UPWIND DIFFERENCE SCHEMES 


In reference 19, Fromm's zero-average-phase-error scheme (ref. 16) for integrat- 
ing the linear convection equation (2-5) was made monotonic by van Leer through the 
inclusion of nonlinear feedback terms (or flux limiters). In this section, we will 
extend van Leer's second-order "monotonicity-preserving" scheme for a single 
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conservation law to a hyperbolic system of conservation laws in a very natural way 
through the use of the normalized Jacobian matrix A* and its associated properties 
described in section 2 . 

A second-order accurate TVD scheme using van Leer's (ref. 19) smoothness monitor 
based on the characteristic flux difference splitting concept can be expressed as 
f o 1 lows : 


U n+1 - U n = -A(A + . . .A. , , n F + A . . .A , , n F) - -7 (I 
3 3 j - ( 1 /z ) J-(l/2) J+(l/2) J+(l/2) 4 


s (e j ))[ (A; +(i / 2) 


AA j+(l/2) A j+(l/2) F ( 1 / 2 ) AA j-(l/2)^ A j-(l/2) F -' 

| (i + s(e J _ l ))[(A+_ (l/2) - xa+. (i/ 2 )) a._ (i/2) f - (a;_ (3/2) 


- AA + . . .)A. , , . F ] + 7 - (I - S(0., ))[(aT 1/ , N + AA. . , . ) A . , , .F 
J"( 3/2)^ J-(3 / 2 ) 4 J+l j+(3/2) l+(3/2)^ J+(3/2) 


( V(x/2) + AA i + (l/2) ) V(l/0 Fl + 4 (I + S(e j ))[( V (l /2) 


+ A A . ,.)A. / . , ,F - (A, , ,.+XA. . ,v)A. . , ,F 

J + (l/2> J + (l/2) J-(l/2) J — ( 1 /2 ) J-(l/2) - 


(4-1) 


Here is a diagonal matrix given by 

S( 6 .) = diag{ s (0 . ) } , 1 = 1 , 2 , . 

J ^ 3 


m 


where 


MV 


r o. 


i 

A i+(i/2) u £ 

- 

A i-(l/2) U fl ) 

t 1 

A j+(l/2) U £ 

+ 

A j— (l/2) U J } 


if lVu/2)“d + l V<i/o u i | 


otherwise 


(4-2) 


The so-called "smoothness monitor" 0j is a quantity that, in some way, measures the 
rate of change of Au^ across a nodal point. It is given by 


A j+(i/2)V A j-(i/2) U £ 


For a more detailed description on this subject, the reader is encouraged to read the 
original papers of van Leer (refs. 19,25,26). Without smoothing (i.e., S£(0j) = 0), 
and for the constant coefficient case, equation (4-1): is equivalent to Fromm's scheme 
and is the simplest upstream-centered scheme of second-order accuracy and may be 
regarded as the average of the central difference scheme of Lax and Wendroff (ref. 3) 
and the purely upwind second-order scheme of Warming and Beam (ref. 17). The flux 
limiter defined by (4-2) is applied to scheme (4-1) exclusively. There are other 
nonlinear smoothing techniques available such as the ones devised by Harten (ref. 18) 
and by Roe (ref. 20). In reference 18 Harten has described a recipe of converting a 
3-point: first-order accurate TVD scheme into a five-point second-order accurate TVD 
scheme by using a modified flux function. In the following, a second-order TVD scheme 
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is described using Harten's (ref. 18) modified flux approach and based on equa- 
tion (2-15). 


First, a modified flux vector is defined for a hyperbolic system of conservation 
laws based on the characteristic flux difference splitting concept. Then, applying 
our first-order scheme to this modified flux vector, we solve 

-M M 

3 U + A 3 F = 0 (4-3) 

t x 

M 

where F is the modified flux vector and its value at nodal point j is given by 

F 1 ? = F. + G./X (4-4) 

1 1 J 

where F is the original flux vector and G is an additional vector remaining to be 


defined. The column vector G at nodal point j 
with its m components given by 

is G. = 
3 

(8l ,j ,g 2,j’ * * * ’ g m,j )T 

8 *,j = ^,j + (i/2) m±n( l g £,j + (i/2)M^,j-(i/ 2 )! ) ’ 

when 

g £,j + (i/ 2 ) g £,j-(i/2) “ 0 

= o , 

when 

g £»j + (l /2) g £,j-(l/2) ~ 0 



(4-5) 


where g^ j + ( x / 2 )^ = 2» . • . , m) are components of the following column vector 

G,,, , . = XsgnA . . .A... , -F/2 - X 2 A.. / , .A,. , , s F/2 (• 

J+(l/2) 6 J+(l/2) J+(l/2) J+(l/2) J+(l/2) 

and 


S «..j + (i/2) Sgn( %,j + (i/2) 
The sgnA in equation (4-6) is given by 


) 


sgnA = Qdiag{ sgn a }Q 


-l 


(4-7) 


(4-8) 


aM M 

and the matrix A in equation (4-3) is related to A which is defined as follows: 


A M = Qdiag{ a £ )Q 1 = SF^/aU 


(4-9) 


A second-order accurate TVD scheme based on characteristic flux difference splitting 
for equation (4-3) is 


IT n+i _ _n . _M -- M. 

3 3 3 ~( 1 / 2 ) J-( 1 / 2 ) 3 +( 1 / 2 ) J+( 1 / 2 ) 


(4-10) 


In equation (4.10) we have assumed that 


-M+ _ ; ± 
A ~ A 


The second-order schemes defined by equations (4-1) and (4-10) are total variation 
nonincreasing which can be proven (at least for the scalar case; see Sweby, ref. 22). 
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5. THIRD-ORDER UPWIND DIFFERENCE SCHEME 


In references 25 and 26, van Leer also investigated third-order schemes for 
numerical convection based on upstream-centered differencing which is regarded as 
higher-order sequels to Godunov's method (ref. 6). Before van Leer, Wesseling 
(ref. 27) derived and tested various high-order upstream schemes, among which are 
Fromm's (ref. 16) scheme, and the third-order scheme. His tests also include several 
central difference schemes, among which are the Lax-Wendroff (ref. 3) scheme and the 
third-order (RBM) scheme derived by Rusanov (ref. 28) and by Burstein and Mirin 
(ref. 29). It is well known that for five-point two time-level upwind schemes the 
highest order which can be achieved is three. It is also found that third-order 
schemes have less dispersion error than the second-order schemes (refs. 26 and 30). 
Thus it is interesting to see how much can be gained in going to third-order accuracy. 

In this section, we describe a third-order upwind (or upstream-centered) scheme 
for integrating hyperbolic systems of conservation laws by applying our characteristic 
flux difference splitting concept to a third-order upstream-centered scheme for the 
linear wave equation (eq. (35) of ref. 25) . The procedure is similar to the one we 
describcid previously for the second-order scheme using van Leer's flux limiter 
(ref. 9). 


+ By using the peculiar properties (eq. (2-24)) of the normalized Jacobian matrices 
A~ and the Jacobian matrices A - of a hyperbolic system of conservation laws, a 
third-order accurate upwind (upstream-centered) difference scheme for integrating 
equation (2-15) can be expressed as follows: 

°r - D " ■ - x<A j-( 1 /oV<i/o F + V(.w‘«./o r) - i (I - s( V> [2J $ + <i/2) 


3AA 


j+(l/2) 


9 +2 

+ a 2 a. 


- (2A + ,, /o , - 3AA + . , . 
J + (l/2) j+ ( l / 2 ) 3 - ( i /^ ) 1— (1/2) 


+ A 2 A +2 , . . )A . , , . F ] - 2 - (I + S(6. ))[(A + , > , 

j-(i/ 2 ) / J — ( 1 / 2 ) 6 3-1 J-(i/z) 


, +2 

a 2 a; 


j-(l/2)^ A j-(l/2) F 


(A j-( 3 / 2 ) A2A j-(3/2) )A j-(3/2) F ^ + 6 (I S(6 j+l ))[(A j + (3/2) 


- A A 


-2 




V F - (A 


A Z A . 


,)A. 


j -+- ( 3 / 2 ) j + (3/2) j+(l/2) j + (l/2) 2 j + (l/2) 


J] 


_2 


+ Y (I + S(0.))[(2A , , , + 3AA. 1 , / n \ + A 2 A. , . . . )A . , .F 

6 J+(l/2) J+(l/2) J+(l/2) 2 J+(l/2) 


(2A. 


■j-(i/2) 


+ 3 AA . 


(1/2) 


+ A A . 


_2 


)A. 


1 “ ( 1 /a ) j-(i/2) 


(5-1) 


where Fs defined as in equation (4-2). Note that all members of equation (5-1) 

are perfect differences; hence, the scheme is conservative. Whether the scheme 
defined by equation (5-1) is TVD or not is yet to be proven. But the flux limiters 
for the present third-order scheme seems to work remarkably well for one-dimensional 
problems; numerical results indicate the scheme is monotonicity-preserving. Disper- 
sive and dissipative properties and the global order of accuracy of this third-order 
scheme will be examined in the following section. In this report, the third-order 
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scheme (5-1) is currently limited to one-dimensional problems only. Present knowledge 
of nonlinear smoothing techniques for third- and higher-order schemes is rather 
limited, and further study of this is warranted (refs. 26 and 31). 

A linear stability analysis, assuming A = constant, yields the following stabil- 
ity bound for scheme (5-1) : 

| A | A < 1 (5-2) 

Equation (5-2) is also the stability bound for the second-order schemes equa- 
tions (4-1) and (4-10). 


6. STABILITY AND ACCURACY OF THE SECOND- AND THIRD-ORDER SCHEMES 


To perform a linear stability analysis, the second-order scheme equation (4-1) 
and the third-order scheme equation (5-1) are applied to the linear wave equa- 
tion (2-5) which is the model equation for the hyperbolic system equation (2-2). 

We rewrite equation (2-5) as 

3 u + a9 u = 0 (6-1) 

t x 

A general expression for explicit two time-level difference schemes which 
approximate equation (6-1) is 


n+i 
u . 

1 


k 2 

E c k u "+k 


( 6 - 2 ) 


With a assumed constant, the amplification factor g of equation (6-1) is defined 
in the customary way as the ratio of the amplitudes of a harmonic wave 
u(x,t) = u(t)exp (ivx) at time t + At and time t. Thus one finds: 

g(0) = exp (-io0) (6-3) 

where a = A a is the Courant number, and 0 = vAx. 


Equation (6-3) has unit modulus and a phase shift per time increment At of 


cp = -avAt = -o0 (6-4) 

e 

The amplification factor g^ of the difference scheme equation (6-2) is found to be: 

k 2 

g h ( 0 ) = ^ c k exp (ike) (6-5) 

k i 


For scheme equation (4-1), it is found that the coefficients c^ are given by 


c = -Aa" l "(l - Aa + )/4 


(6-6a) 
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c i = Aa + 4- Aa(l - A|a|)/4 

(6-6b) 

c 0 = 1 - 3A | a | /4 - A 2 a 2 /4 

(6-6c) 

c +i = -Aa - Aa(l - A|a|)/4 

(6-6d) 

c, = Aa (1 4- Aa )/4 
+2 

(6-6 e) 


With these coefficients and equation (6-2), it is easy to find that the Taylor exapn- 
sion of the scheme is 

Au 1 ! = u n+1 - u° = -AaAx(3 u) . + —• (AaAx) 2 (3 u) . 

3 3 3 x 'j 2 xx 'j 

4- (AaAx 3 / 1 2 - A 2 Ax 2 a j a | ) ( 3 u) . (6-7) 

XXX j 

which is of second-order accuracy both in time and space. 

The amplification factor for the second-order scheme is 
g^(6) = Re g (0) + i Im g^(0) = 1 - [ cr | ( 1 - cos 8)[1 4- |a| - (1 - |a|)cos 0]/2 

- io sin 0[3 - [ a | - (1 - |o|)cos 0]/2 (6-8) 

The phase shift <J> per time step of the finite-difference scheme can be computed by 

the formula 

<P = sin - 1 { Im g h ( 0 ) / I g h I > (6-9) 

The relative phase error of scheme (4-1) is given by the ratio of (6-9) and (6-4). 

Dissipative (amplitude distortion) and dispersive (phase shift) errors can be 
conveniently represented as polar plots of | gp ( 0 ) ! and <j>/<p e , as a function of 0. 

The same procedure can be carried out for the third-order scheme (5-1). 


The coefficients c^ for (5-1) are given by 

c_ 2 = -Aa + (1 - A 2 (a + ) 2 )/6 (6-10a) 

c = A | a | ( 1 - A 2 a 2 ) /3 + Aa(l + Aa) (4 - Aa)/6 (6-10b) 

c 0 = (1 - A | a | /2 ) (1 - A 2 a 2 ) (6-10c) 

c = A | a | ( 1 - A 2 a 2 )/3 - Aa(l - Aa) (4 4- Aa)/6 (6-10d) 

c +2 = Aa~ (1 - A 2 (a _ ) 2 ) /6 (6-10e) 

The Taylor expansion for equation (6-2) with (6-10) is 
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Au? = u n+1 - u 1 ? = -AaAx(3 u) . + — (AaAx) 2 (8 u) . - — (AaAx) 3 (3 u) . 

J j 3 x '3 2 xx '3 6 xxx '3 

- ± A| a| (1 - A | a | /2 - A 2 a| a| )Ax 4 O^u) . (6-11) 

The accuracy of this scheme is third-order in time and space. 

The amplification factor for the third-order scheme is given by 

g h (6) = 1 - | ct | ( 1 — cos 6 ) [ 1 + 3 1 ct | - a 2 - (1 - |a| 2 )cos 0]/3 

+ ia sin 0 [ ( 1 - a 2 )cos 0 - (4 - a 2 )]/3 (6-12) 

The dissipation and dispersion of (6-12) can be similarly computed. 

Figure 2 shows plots of the amplification factor modulus and phase error of the 
second- and third-order schemes for several values of the Courant number o in the 
stable range. 

In terms of dissipation, the third-order scheme is slightly less dissipative than 
the second-order scheme. But also, in terms of disperson, the third-order scheme has 
a larger dispersion. 

In order to check the global order of accuracy of the second-order and the third- 
order schemes, we next describe some numerical experiments with periodic scalar con- 


servation law 

3 u + 3 (u 2 /2) =0 -it < x < 7T , t > 0 (6-13a) 

t X 

with initial data 

u(x,0) = 2 + sin x ; -tt < x < tt (6-13b) 

and the periodic boundary condition 

u( — TT , t ) = U ( TT , t ) , t > 0 (6-13c) 

It is easy to see that the solution to equation (6-13) is smooth up to t = 1 , at 
which time a shock is formed (see ref. 32 for a more detailed description). The 


relative L 2 -errors are computed at the same physical time for the second-order 
schemes (4-1) and (4-10) and the third-order scheme (5-1). 

In table 1 we present a mesh refinement chart for the schemes defined by (4-1), 
(4-10), and (5-1). For comparison purposes, we also list the results obtained by 
using the central difference scheme of Lax and Wendroff (LW) (ref. 3), and the second- 
order purely upwind difference scheme of Warming and Beam (WB) (ref. 17). The second- 
order scheme (4-1) without flux limiter, denoted by (II), is identical to the zero- 
average phase-error scheme of Fromm (F) (ref. 16). (Il-a) denotes the second-order 
TVD scheme with van Leer's flux limiter equation (4-1). The scheme using Harten's 
(ref. 18) modified flux defined by equation (4-10) is denoted by (Il-b). The third- 
order schemes defined by equation (5-1) without and with flux limiters are represented 
by (III) and (Ill-a) , respectively. 
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The values Ax and At in table 1 are defined by 

Ax = 2 tt /M , At = 0. 95Ax/3 

The relative L 2 -error in this chart correspond to results at t = 0.497, at which 
the solution is still smooth. The relative L 2 -error used in table 1 is defined by 



where uj is the numerical solution and Vj is the exact solution obtained by 
Newton-Raphson iteration. 


From the data of the numerical experiments in table 1, we observe that the 
second- and third-order schemes without flux limiters (II and III) are fully second- 
and third-order accurate; a refinement by a factor of 2 reduces the error by a factor 
of approximately 4=2 and 8=2, respectively. With flux limiters, the second- 
and third-order schemes have about the same quality of results; the global order of 
accuracy for both schemes is around 3/2. In general, the upstream-centered difference 
schemes (Il-a and Ill-a) demonstrate more accurate than the central difference schemes 
(LW and Il-b) . 


7. THE EULER EQUATIONS OF INVISCID GAS DYNAMICS 


In this section, we apply the above development to the Euler equations of invis- 
cid gas dynamics in one and two space dimensions. The one-dimensional Euler equations 
can be written in conservation-law form as 


3 U + 3 F (U) = 0 (7-1) 

t x 

T 

where U = (p,pu,e) is the state vector of conservative variables and 
F = (pu,pu 2 + p,(e + p)u)^ is the flux vector. The so-called primitive variables of 
equation (7-1) are the fluid density p, the velocity u, and the pressure p. For a 
perfect gas, the total energy per unit volume, e, is related to the other variables 
by e = p/ (y - 1) + 0 .5 pu“, where y is the ratio of the specific heats. Equa- 
tion (7-1) can be expressed in nonconservative quasi-linear form as 


3 U + A3 U = 0 
t x 

where A is the Jacobian matrix 3F/3U 


0 1 0 
(y - 3)u 2 /2 (3 - y)u (y - 1) 

(y - l)u 3 - yeu/p ye/p - 3(y - l)u 2 /2 yu 

. 

The eigenvalues of A are 


(7-2) 


(7-3) 


= u 


= u + c 


and 


= u 


(7-4) 
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where c = /yp/p is the local speed of sound. The coefficient matrix A can be 
diagonalized by the similarity transformation matrix Q = MT with M and T given by 




■ 1 

0 

3U/3U 

= 

u 

P 



_u 2 12 

pu 1/ 


"l 

p /2c 

-p /2 c" 

T = 

0 

1/2 

1/2 


.0 

pc/2 

-pc/2_ 


0 

0 


(7-5) 


(7-6) 


T -l 

where U = (p,u,p) . M and T are the inverse matrices of M and T and can be 
easily obtained. 

We have 


A = T 1 M~ 1 AMT 


u 0 0 

0 u + c 0 

0 0 u - c 


The matrix A is split into two parts as 


.+ 


+ 


A = A' + A = MT ( A ' + A )T -1 M -1 


Similarly, we have 




+ 


■i«-i 


A = A + A = MT (A + A )T M 


As indicated in equation (2-22) we can also evaluate ^•_( 1 / 2 ) as 

v + 

V j-(l/2) 


A 1 , , , = MTA~T 1 M~ 1 


where M and T are given by 


M = 


1 0 0 
up 0 

|_(u) 2 /2 pG 1 / (y - 1)J 


and 


T = 


1 p/ (2c) -p/ (2c) 


0 1/2 1/2 

LO pc/2 -pc/2 


(7-7) 


(7-8) 


(7-9) 


(7-10) 


(7-11) 


(7-12) 
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The overbar quantities are evaluated by a two-point average. 


For example, (u)j_(i/ 2 ) is defined as (hj _ 1 + uj)/2. For the Euler equations in 
multiple space dimensions, similar procedures can be applied; below we briefly 
describe some of the highlights. Details of the formulation for the three-dimensional 
Euler equations in general curvilinear coordinates can be found in reference 14. 


The unsteady two-dimensional Euler equations in conservation-law form are given 
by 

3 t U + 9 X F + 9 y G = 0 (7-13) 

where U = (p ,pu,pv,e) T , F = (pu,pu 2 + p,puv,(e + p)u) T , and 
G = (pv,pvu,pv 2 + p,(e + p) v) . 


In the above equation p is the density, u and v the velocity components 
corresponding to the x and y coordinates, e the total internal energy per unit 
volume, p the pressure, and y the specific heat ratio. 


The equation of state is 

P = (Y - 1) |e - j p (u 2 + v 2 ) 


>] 


Applying the characteristic flux difference splitting method to equation (7-13) we 
have 


9 t U + (A + + A ) 3 x F + (B + + B ) 9 y G = 0 


where 


A - = MT A - T -1 M~ 1 and B“ = MT A i T -1 M' 


L 


XXX 


Matrices M, T^ 1 , and T 1 are given by 


- 1 M -1 


M" 


1 

-u/p 

-v/p 


0 

1/P 

0 


y y y 

o 

o 

i/p 


o 

o 

o 


(y - 1 ) (u 2 + v 2 )/2 -(y - l)u -(y - l)v (y - 1 ) 




0 0 


P 0 


0 1 


P o 


-1/ (c) 2 

l/s 

0 

-1 /c . 


(7-14) 


(7-15) 


(7-16) 


(7-17) 
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and 


10 0 1 / (c ) 2 "1 



0 10 0 

0 0 p 1 /c 


,0 0 p -1/3 J 


(7-18) 


A second-order scheme using van Leer's (ref. 19) flux limiter can be constructed 
as in the one-dimensional case, applying the same flux difference components for F 
to the flux vector G and simply replacing the A's by the matrices B's. Numerical 
experiments for two-dimensional Euler equations using this scheme have been reported 
in reference 12. 

Similarly, we can apply Harten's (ref. 8) modified flux to solve 

3 U + (A + + A") 3 F M + (B + + B _ )3 G M = 0 (7-19) 

t x y 

M 

The modified flux in the x direction F is defined in the same manner as in 
Section 4. The modified flux in the y direction G 1 can be similarly defined. 

An explicit second-order TVD scheme based on equation (7-19) can be expressed as 


JM 


U n+ ^ - U n = -A[A + , . (F M - F M ) + A~ . . . (F\ , 

J,k 3 , k j -(1/2) ,k j,k j-i,k y j+(i/2),k^ 3+1, k 


M 

O 1 


A[B I,k-(i/0 (G j,k 


M M 

G . . ) + B . , , . . . (G . , , 

j,k+(i/2) v 3 ,k+i 


G^ k )] (7-20) 


For implicit schemes, Strang-type dimensional splitting can be applied to construct 
robust schemes as described in reference 14. 


For flow problems involving complex geometry, the above development for the 
Euler equations in Cartesian coordinates can be easily extended to general curvilinear 
coordinates. For a multidimensional case, the third-order scheme with the flux 
limiter defined by equation (4-2) is not sufficient to make the scheme stable, further 
study on the flux limiter (nonlinear filtering techniques) for higher order "mono- 
tonic" upwind schemes is needed (ref. 31). 


8. RESULTS AND DISCUSSIONS 


In this section, we present the results of some numerical experiments using the 
second- and third-order upwind difference schemes which we have described previously. 
Single conservation laws and systems are considered. For single equations, we con- 
sider the inviscid Burgers' equation, and, for systems of equations, the Euler equa- 
tions of inviscid gas dynamics in one and two space dimensions are considered. For 
third-order schemes, only results for one-dimensional problems are presented. For all 
the following numerical calculations, uniform computational mesh systems are employed. 


18 



(a) Inviscid Burgers' Equation 


The first model problem we consider is the inviscid Burgers' equation. 

3 u + 3 x f(u) = ® ’ f( u ) = u 2 /2 

with initial conditions 


u(x < 0) = u^ , and u(x > 0) = u^ 


The ul and up values can be chosen to represent different situations which may 
occur in transonic aerodynamics. We consider the following four cases: 


(a) transonic moving shock: u^ = 2 , u R = -1 

(b) transonic stationary shock: = 1 , u = -1 

(c) transonic expansion: u = -1, u = 1 

L K 

(d) fully supersonic: u = 0.5, u^ = 1 

The following initial values were prescribed 


u = u^ for j < J 


u.=u„ for i > J + 1 

J R 

For cases (a-c) , J is equal to 25 and for case (d) J is 16. The total number of 
mesh points is 50 and the mesh size is Ax = 0.01. Figures 3(a)-3(c) display the 
solutions using scheme (2-5a,b) for cases (a-c). Transonic moving and stationary 
shocks are captured within two zones and the transonic expansion shock is excluded 


Figures 4(a) and 4(b) show results for case (c) using the second- and third- 
order accurate schemes. 


For case (d), we solve the case as reported by van Leer in reference 19 for com- 
parison purposes. Figures 5(a)-5(d) show the second- and third-order results without 
and with nonlinear smoothing, respectively. 

In order to further test the present method we also include the solution to 
equation (6-13) at two instants; one at time t = 0.497, at which the solution is 
still smooth and another at time t = 0.995, at which time a shock is formed. 

Figures 6(a), 6(b), 7(a), and 7(b) show the results using the second-and third-order 
schemes (II, Ila, III, and Ilia) with 20 grid points (M = 20 in table 1). 


(b) Ideal Shock Tube Problem (One-Dimensional Euler Equations) 

We consider the ideal shock tube problem as our second model problem for the 
one-dimensional Euler equations. The case we consider here is the one reported by 
Sod (ref. 27). The initial conditions defined in the problem are 
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P L = 1 , P R = 1 , ^ = 0 , 0<x<0.5 

p R = 0.1 , p R = 0.125 , u^ = 0 , 0.5 < x < 1 

The space-step size Ax used is 0.01. The time step At is calculated at each time 
level subject to the condition that the maximum Courant number is 0.9, and the output 
is at time t = 0.24. The velocity, density, pressure, internal energy, and entropy 
profiles are presented for all calculations. The numerical solution (box with con- 
necting line, every other point) and the exact solution (solid line) are plotted. 

Figures 8(a)-8(e) indicate the results using scheme equation (4-1) without flux 
limiters (II). Figures 9(a)-9(e) show the results using scheme (4-1) with flux 
limiters (Ila) . These results have been reported in reference 14 and are comparable 
to those reported by van Leer using his MUSCL scheme (ref. 34). The MUSCL scheme is 
a Lagrangean scheme, and the Lagrangean results are remapped with least-squares 
accuracy onto the desired Euler grid in a separate step. The present scheme (Ila) 
can be viewed as the corresponding Eulerian scheme for MUSCL scheme. 

Figures 10(a)-10(e) display the results using equation (4-10) with Harten's 
(ref. 18) modified flux (lib) . These results invite comparison with those obtained 
using Harten's high resolution scheme (ref. 18). 

Figures ll(a)-ll(e) show the third-order results using scheme equation (5-1) with 
van Leer's (ref. 19) flux limiters (Ilia). 

It is observed that the shock profile is sharpened compared with second-order 
results. For the contact surface profiles, the one using van Leer's (ref. 19) flux 
limiter is slightly less smeared than that using Harten's (ref. 23) modified flux. 

To appreciate the effect of flux limiters, the third-order results without flux 
limiters (III) are also shown in figures 12(a)-12(e). The flux limiters defined by 
equation (4-2) seem to work very well in this one-dimensional problem. 


(c) Quasi-One-Dimensional Nozzle Problem 

The nozzle geometry we consider here is a divergent nozzle with area distribution 

A(x) = 1.398 + 0.347 tanh(0.8 x - 4) , 0 < x < 10 

We consider the case of supersonic inflow and subsonic outflow with a stationary 
shock. For all explicit schemes, we use CFL =0.9. Uniform mesh distribution with 
mesh size Ax = 0.2 (50 grid points) is used. The converged solution is obtained in 
about 800 iterations (with maximum residual less than 10~ 6 ). The density and mass- 
flux distributions are plotted. 

Figures 13(a), 13(b), 14(a), and 14(b) show results using the first-order 
explicit scheme and averages defined by equation (2-22) without and with transonic 
shock operators, respectively. We can see that the spike which appears near the 
post-shock region in figure 13 has been eliminated in figure 14. 

In figures 15(a) and 15(b), we show the corresponding results using the first- 
order implicit scheme with averages defined by equation (2-23) and CFL = 25. 
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Both figures 14 and 15 display a clean shock profile which occupies one mesh 
interval. 

In figures 16 and 17, we show results using the second-order explicit schemes 
(Ila) and (lib), respectively. Comparing with the first-order results, the second- 
order effect can be seen from the mass-flux profile near the post-shock region. 

For steady shock calculations, we found that there is no need to apply transonic 
shock operators if one uses the averages defined by equation (2-23) (here eqs. (7-10) 
to (7-12)). This type of average, in general, captures steady-state shocks within 
one transition cell as also observed by other approaches (refs. 11 and 8). 

In figures 18(b) and 18(b), we present the results using the third-order scheme 
equation (5-1). 


(d) Plane Shock Reflection (Two-Dimensional Euler Equations) 

In order to test the current characteristic flux difference splitting method in 
a multidimensional problem, we consider the plane shock reflection problem. Here a 
simple regular reflection and a Mach reflection are considered. Preliminary results 
using second-order schemes only are presented. For the case of regular reflection 
shock, we solve the case reported by Yee et al. [24]. The incident shock angle is 
29° and the free stream Mach number is 2.9. The computational domain is 0 < x < 4, 
and -0 < y < 1, with uniform grid size of 61*22. 

Figure 19(a) and (b) show the density contour and the pressure distribution at 
y = 0.5 using scheme Ila. Figures 20(a) and (b) show the same results using 
scheme lib. For the case of simple Mach reflection, the incident shock angle is 42° 
and the free-stream Mach number is 4.23. The computational grid is 34 by 48 with 
Ax = Ay = 0.1. Density and pressure contours are shown to demonstrate the resolution 
and the physical structures as well. Figures 21(a) and (b) show the results using 
scheme equation (7-19) with Harten's (ref. 23) modified flux approach. The results 
using the implicit version of scheme (4-1) were reported in reference 12. Further 
tests and results for more complex flow problems using generalized coordinate systems 
will be reported in a future study. 


9. CONCLUSIONS 


In this study, we have described second- and third-order two time-level and five- 
point upwind difference schemes for hyperbolic systems of conservation laws and their 
application to the Euler equations of inviscid gas dynamics. Entropy-satisfying 
transonic shock transition operators are introduced at the sonic point in accord with 
the principle of upwind differencing which enable the schemes to exclude nonphysical 
expansion shocks. Stability and accuracy of the second- and third-order schemes are 
examined. Nonlinear smoothing techniques devised by van Leer (ref. 19) and by Harten 
(ref. 23) are adapted and generalized to the present framework of characteristic flux 
difference splitting to make the schemes total variation diminishing (at least for a 
scalar equation). Numerical experiments also justified this claim. The method offers 
some simplicity in terms of the mathematical formulations and its accordance with the 
idea of preserving the physical domain of dependence, hence is adequate for the 
numerical solution of partial differential equations governing convection (ref. 35). 
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For the second-order schemes, our results are comparable with those obtained by 
using van Leer's (ref. 19) and Harten's (ref. 23) high resolution schemes. The third 
order scheme, although only tested for one-dimensional problems, displays slightly 
superior resolution compared with second-order results. Whether or not the present 
third-order scheme is TVD is yet to be proven; but the results seem to indicate that 
the scheme is monotonicity-preserving. Further study on higher-order flux limiters 
for integrating convective-type equations is suggested. 
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TABLE 1.- COMPARISON OF L 2 -ERROR FOR SECOND- AND THIRD-ORDER SCHEMES 


Ax 

LW 

WB 

F (or II) 

Il-b 

Il-a 

III 

Ill-a 

2tt 

20 

2 TT 
40 

2 TT 
80 

2 TT 
160 

0.449E-01 
. 128E-01 
. 338E-02 
. 863E-03 

0.380E-01 
. 11 IE-01 
.294E-02 
. 754E-03 

0 .2 17E-01 
.488E-01 
. 104E-02 
. 230E-03 

0.295E-01 
. 943E-02 
. 308E-02 
. 973E-03 

0 . 204E-01 
. 620E-02 
. 177E-02 
. 503E-03 

0.144E-01 
•241E-02 
. 338E-03 
. 442E-04 

0. 166E-01 
.543E-02 
. 174E-02 
.534E-03 




Figure 1.- Transonic shock transitions. (a) Transonic shock; (b) transonic 

expansion. 
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Figure 5.- Concluded. 
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Figure 6.- Solution of the Burgers equation. (a) Smooth solution (II); (b) smooth solution (Ila) ; 

(c) shock solution (II) ; (d) shock solution (Ila) . 



Figure 7.- Solution of the Burgers equation. (a) Smooth solution (III); (b) smooth solution (Ilia); 

(c) shock solution (III) ; (d) shock solution (Ilia) . 
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Figure 9.- Ideal shock tube solution (second-order, Ila) . (a) Density 

(b) velocity; (c) pressure; (d) internal energy; (e) entropy. 
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Figure 10.- Ideal shock tube solution (second-order, lib) . (a) Density; 
(b) velocity; (c) pressure; (d) internal energy; (e) entropy. 
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Figure 14. Divergent nozzle (explicit first-order, with shock operator). 

(a) Density; (b) mass flux. 
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Figure 17.- Divergent nozzle (explicit second-order, with shock operator). 

(a) Density; (b) mass flux. 
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Figure 19.- Regular shock reflection (scheme Ila) . (a) Density contours 

(b) pressure distribution at y = 0.5. 
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Figure 20.- Regular shock reflection (scheme lib) . (a) Density contours; 

(b) pressure distribution at y = 0.5. 
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Figure 21.- Simple Mach reflection (scheme Ila) . (a) Pressure contours; 

(b) density contours. 
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